The Hartree-Fock ground state of the three-dimensional electron gas 
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In 1962, Overhauser showed that within Hartree-Fock (HF) the electron gas is unstable to a spin 
density wave (SDW) instability. Determining the true HF ground state has remained a challenge. 
Using numerical calculations for finite systems and analytic techniques, we study the HF ground 
state of the 3D electron gas. At high density, we find broken spin symmetry states with a nearly 
constant charge density. Unlike previously discussed spin wave states, the observed wave vector of 
the SDW is smaller than 2kF- The broken-symmetry state originates from pairing instabilities at 
the Fermi surface, a model for which is proposed. 
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The three-dimensional electron gas is one of the basic 
models of many-body physics, and has been investigated 
for over 70 years P, % H, S, S IS H- As the simplest 
model system representing an itinerant metal, it consists 
of interacting electrons in a uniform neutralizing charge 
background, described by the Hamiltonian: 
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where the sums are over particle indices. Its properties 
are routinely used in density functional theory, e.g., in 
local density approximations, as a reference state in cal- 
culations of electronic structure of real materials • 

The simplest approach to an interacting many-fcrmion 
system such as the electron gas (jellium) is the mean-field 
Hartree-Fock (HF) method, which finds the Slater deter- 
minant wave function minimizing the variational energy. 
In unpolarized jellium, the "conventional" solution is a 
paramagnetic state with spin symmetry, the restricted 
HF (rHF) solution in quantum chemistry. 

The rHF solution, however, is not the exact HF ground 
state of jellium. In 1962, Overhauser Q proved that the 
rHF solution is unstable with respect to spin and charge 
fluctuations at any density. The global minimum en- 
ergy state within HF is a spontaneously broken symme- 
try state. The properties of this global ground state have 
remained unknown Q. This is surprising, given the fun- 
damental importance of both the electron gas and the HF 
approach. The correlation energy of the homogeneous 
electron gas is a commonly used fundamental concept, 
but its definition is in terms of the HF energy of the 
electron gas. 

In this paper, we numerically find the HF ground state 
for finite systems. Our motivation, aside from solving 
this mathematical puzzle, was to understand the mech- 
anism for the broken symmetry state. Further, it was 
hoped that the solution would suggest candidate ground 
states, for jellium or for other systems, that can then 
be explored by accurate many-body approaches such as 
quantum Monte Carlo [1, H]. We focus on high and 



medium densities. An analytic approach is used to aug- 
ment and extend the results to the thermodynamic limit. 
We find that a different pairing instability characterizes 
the high-density ground state. 

We consider N (iV| = iV| = N/2) electrons in a cu- 
bic supercell of volume Vl = . The density is spec- 
ified by the average distance between electrons: = 
{i^/ AttNY^^ / as- We write a Slater determinant as 
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with = X)k ''i C*^)!^) where |k) is a plane-wave basis 
function. In the rHF solution, any |k) with A: = |k| < kp 
is fully occupied, while all others are empty. Our basis 
contains all plane- waves with k < fccut ■ 

To find the global ground state, i.e., the unrestricted 
HF (uHF) solution, wc use an iterative projection 
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where iJHF(^'™'') is the HF Hamiltonian, i.e., the mean- 
field approximation of Eq. ([1]). The wave function re- 
mains a single Slater determinant in the projection 
If T is sufficiently small, the energy will decrease in each 
step and the projection will converge to a HF solution as 
m oo. To ensure that the solution is not a local mini- 
mum, wc often start from multiple random initial states 
$'"^) and verify that the same final state is reached. 

The smallness of the energy scale relevant to the sym- 
metry breaking at high density presents a difficulty; 
finite-size effects can easily be larger than differences in 
energy of different phases. As a consequence, the stable 
structures vary wildly with N. There are subtle com- 
mcnsuration effects in both r-space (Wigner crystal) and 
k-space (spin waves). For example, in open-shell systems 
there is always a uHF solution, since a broken symmetry 
state can be formed in a partially filled shell to lower the 
exchange energy, with no cost to the kinetic or Hartree 
energy (see below). In closed-shell systems, on the other 
hand, there seems to be a critical value of r^{N), below 
which no uHF state exists for a given value of N under 
periodic boundary condition (PEC). 



To break the shell structure, we impose twisted bound- 
ary conditions Q on the orbitals (and hence on the wave- 
function): (j>{r + La) = e*^'^*°(/)(r) where a = x,y,z. 
This applies a shift of kg = 2TT9/L to the plane- wave ba- 
sis. 9 = corresponds to PBC, the F-point for solids. 
For generic twist angles 6, the rHF solution is non- 
degenerate. 

Let us write the Hamiltonian in second quantized form, 
omitting an overall constant: 

(T.k A ^ 

where c^^ ^ and Ck.o- are creation and annihilation op- 
erators. In the second term, A denotes the variables 
{k, k', Q, (T, a'}, Q is a reciprocal lattice vector, the ' on 
the summation indicates Q 0, and 

^(A) = cLQ,^4'-i-Q,<T'Ck',^'Ck,<T- (5) 

The HF Hamiltonian 77hf(«'^"^) needed in Eq. © is 
the same as H, but with V{A) replaced by the linearized 
form Vhf(A) = ^(A) - (u(A)}/2, where 

*(A)ee 2[(Ek"4",.'Ck"-Q,.')4',k+Q (6) 

^(4'+Q,<T'Ck+Q,<T')5<T,^' ] cl aCk\a- 

The expectation (..) is with respect to |$("'^). The vari- 
ational energy is £'t,(<i>) = (<l>|if|$)/($|$). We also 
compute a "growth estimator" of the energy, Eg = 
_ln[($(™+i)|$(™+i)}/($(™)|$(™))]/2T, in the projec- 
tion. At convergence, = Eg, 
which means that |$(™)) is an eigenfunction of i^HF- 
Hence, the projection gives a true solution of the HF 
Hamiltonian, not just a Slater determinant with a varia- 
tional energy lower than the rHF value. 

Calculations were carried out with a set of random 9 
and the results averaged and errors estimated. For each 
N, a fixed set of kg-points were used at different values 
of Ts to correlate the runs. At larger r^, less statistical 
accuracy is needed and fewer kg-points were used. Typi- 
cally the plane- wave basis cutoff was set to fccut 
i.e., a kinetic energy cutoff of 5-10 £'f. The resulting ba- 
sis set error is negligible for all but the largest rs (= 7). 
Fast Fourier transform (FFT) techniques were used to 
speed up each step in the projection [8|, and the orbitals 
re-orthonormalized as necessary. 

The top panel of Fig. [1] shows the energy difference SE 
between the uHF ground state and the rHF state. As is 
reduced, the magnitude of SE decreases before becoming 
nearly flat at high densities (r^ < 3) for each value of N. 
The energy lowering remains finite across r^, showing 
that the broken symmetry uHF solution exists for all 
densities, consistent with Overhauser's proof Q. At low 
Tg, SE is roughly — 10~^Ry in the large-sized systems. 
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FIG. 1: (color online) Energy differences (in Ry) between the 
HF ground state and the rHF state, (a) The energy lowering 
per electron, SE = {E — Eihf)/N, vs. Ts for different values of 
N (b) The three components of SE, kinetic (K), Hartree 
(Vh), and exchange (Vex), for A'' = 54 (solid lines) and = 
66 (dashed lines). The inset is an enlargement of the low 
Vs. Error bars are estimated from the results for various hg- 
points. 



Calculations for larger systems are needed to clarify its 
behavior at small r^. The energy differences are very 
small compared to the rHF energy which, at iV ^ oo, is 
ErHp/N = (2.21/r2 - 0.916/rs)Ry: the relative energy 
reduction vanishes as Ts ^ 0. From the bottom panel, 
we see that the lower total energy of the uHF ground 
state is achieved by reducing the exchange energy at the 
cost of increased kinetic energy (exciting electrons into 
k) states with k > kp). The Hartree energy remains 
unchanged up to rs 3. 

We find that the momentum distribution in the uHF 
solution is spin-independent, i.e., n^ Ck) = rix(^) 
broken symmetry in the broken symmetry ground state. 
This is illustrated in Fig. [5] for = 3, but holds in all 
our calculations when fully converged, at all and N . 

Figure [2] also shows that, as decreases, the occu- 
pancy of fc > kp states becomes less pronounced, and 
significant modifications to the Fermi sphere become in- 
creasingly confined to the immediate vicinity of the Fermi 
surface (FS). At small Vs, such modifications tend to be 
principally single pairing states. An example is seen at 
= 1 in = 54: two plane- wave vectors are involved, 
such that a pair of rHF orbitals = 10^-) = |k) become 

= Ck|k) ± Ck^k') (7) 

where k < kp and fc' > kp, and |ckp -I- |ck'P = 1. Such a 
pairing state by itself forms a linear spin-density wave 
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FIG. 2: (color online) Momentum distribution n(k) at dif- 
ferent Ts VS. k/kp- The main graph has A'^ = 54, with 
e = (-0.368,0.172,-0.364). The arrow indicates kp. The 
inset shows [n(k) — nrHF(k)] at = 1 for A'' = 246, with 
9 = (-0.494,-0.425,0.144). Note the primary spike at kp, 
and the paired smaller spikes with one on each side of kp. 
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FIG. 3: (color online) Left: Peak values and locations of the 
Fourier transforms of spin and charge densities for different 
values of in = 54. Errors are estimated from the values 
at different k^-points. Note the logarithmic scale in S{q). 
Right: Contour plot of the spin density a{x,y,z) for a slice 
parallel to the x-y plane. The system has A*' — 512 and — 1. 
The q- vector [263] has the leading So-(q) value, followed by 
[333], and then [4l5] with Scr four times smaller. 



(SDW), with constant charge density and unchanged 
Hartree energy. As the inset in Fig. [2] shows, additional 
pairing states form in the uHF solution involving wave 
vectors further from the FS. Although their amplitudes 
become very small as is reduced, these are important 
to the true uHF ground state, as we discuss later. 

Real-space properties are examined in Fig. [31 The 
charge density is p{r) ~ ri|(r)+n|(r) and the spin density 
is a{r) = 'T.t(r) — ni{r). We measure their Fourier trans- 
forms, e.g., Sp{q) = I / /9(r)e'i-'"dr| ViV. At r, = 7, the 
A'^ = 54 system is an antifcrromagnetic bcc Wigner crys- 
tal. As Ts decreases, electrons become less localized and 
less particle-like. Fluctuation in the charge density be- 
comes much smaller, as indicated by the rapid reduction 
in Sp{q). This is consistent with the vanishing Hartree 
energy in Fig. [TJ Although Sa{ci) also decreases with 
rs, it is much larger and spin symmetry remains bro- 
ken. At rs = 4, the SDW no longer has a bcc structure, 
and its symmetry between a;, y, and z is broken. At 
small rs, the electrons are highly delocalized and wave- 
like. Charge variations are effectively compensated for 
by spatial "double occupancy" of t and | electrons, as in 
the pairing state discussed in Fig. [21 

We measure the characteristic wave vector of the spin 
or charge density wave by g = |q| , where q is the peak po- 
sition of S{q). A wave vector of ^ 2kp seems to have 
always been assumed in previous investigations of the 
SDW states [3, 0|- However, we find the maximum spin 
ordering is at smaller wave vectors, as shown in Fig. [3l 
for N = 54. Consistent results are seen for larger TV, e.g. 
at rs = 2, q^/kp = 1.2(2), 1.4(2), 1.5(3), and 1.5(2) for 
N = 66, 128, 246, and 528, respectively. 

We next prove analytically that an SDW instability 
whose wave vector is < 2kp indeed exists. Let us con- 



sider a system of large but finite iV, with kg = 0. From 
the rHF reference state, we create a broken symmetry 
state with two pairing orbitals as in Eq. ([7]), using two 
points at the FS, k (a highest occupied state) and k' (a 
lowest unoccupied state), as illustrated in Fig. [H The 
energy cost consists of kinetic and exchange terms : 
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where kp = l/(ars), with a — (4/97r)^/'^. Our choice 
of k and k' gives: Ak = |k'| - |k| - {l/2lp){2Tr/ L) = 
l/{2l'^ars), where Ip is defined by kp = Ip {2tt/L), i.e., 
Ip = (3/87r)^/^A^^/^. For fixed rs, the exchange term 
dominates if N is sufficiently large. Choosing Ck and Ck' 
to be real and of 0(1), we can write Eq. ([9]) as: 
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We now create a "satellite" pairing state in the vicinity 
of kandk', i.e., withk'-k = s'-sand |s-k| ~ 0{2n/L). 
We choose the excitation amplitude to be of the particu- 
lar form: Cs' ^ l/(ln/F)^^'' {5 > 0), and thereby Cs ~ 1. 
Using Eq. ([9]) and noting that the difference in the mag- 
nitude of the wave vectors. As = |s'| — |s|, is less than 
the size of circles in Fig. [4l we obtain an upper bound to 
the energy cost for creating the {s, s'}-pairing state 
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The decrease in exchange energy because of "constructive 
interference" between the two parallel pairs is: 
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FIG. 4: (color online) Cartoon of pairing state with < 
2kF- The primary pairing state {k, k'} is at the FS. The 
satellite pairing state {s,s'} can be in the shaded areas, 
within 0{2n/L) of k and k'. Dashed arrow lines illustrate 
Overhauser[3, [3 pairing at 2k f- 
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For sufficiently large Ip, i.e., a large enough system size, 
AV^^^I can always be made larger than the energy costs 
in Eqs. ([T0|) and pT|) . Hence this is an SDW state with 
lower energy than the rHF state. 

The wave vector of the constructed SDW is Qo- = k'— k. 
As Fig. [4] shows, does not need to be 2fci?. The angle 
between k' and q^, 9, can range from [q^ = 2kp) to 
7r/2 {q„ = 0). As 6 increases, more {k, k'} pairing states 
become available on the FS, while the number of possi- 
ble satellite pairs, i.e., the volume of the shaded areas, 
decreases. The optimal choice would be in between. In 
fact, as a crude estimate, the number of {k, k'} pairs is 
oc 27rfc|,sin0, and the number of {s,s'} pairs for each is 
oc (tt — 6). Maximizing their product gives q^ ^ 1.52kp, 
which is consistent with our data. 

In previous approaches 0, 0], pairing is constructed 
from orbitals directly across the FS. The energy lowering 
is driven by interference between such pairs (dashed ar- 
rows in Fig. |4]), which in our model belong to primary 
pairing states. Our approach differs by including the 
satellite pairing states. The interference between the pri- 
mary and satellite states is what makes a general q^ pos- 
sible. This model is supported by the exact numerical 
data in Fig. [2] Clearly, the true ground state goes be- 
yond this model: the energy will be further lowered by 
having more {s, s'} pairs and multiple {k, k'} states, etc. 

We have shown that the uHF states at intermediate 
and high densities have nearly constant charge density. 
They are wave-like, and arise from pairing between states 
on the FS separated by distance g^. The momentum 
distribution is spin-independent. Its deviation from the 
Fermi sea is increasingly confined to the vicinity of the 
FS as rs 0. The SDW wave vector is q„ ~ 1.5{2)kp, 
and its structures are determined primarily by the short- 
range exchange potential In the rHF solution, the 
size of the exchange hole is = 2.3Ars 0. A linear SDW 



with a wave length (i.e., characteristic like-spin separa- 
tion) of has wave vector q = lAOkp. As 0, the 
numerical results are sensitive to the detailed topology 
of the FS in the finite-size systems. The outcome of the 
SDW structure can vary greatly, as it is a delicate balance 
to optimize pairing among a small number of plane- wave 
states at the FS that can participate. Our results indi- 
cate that, at high densities, the HF ground state tends 
to further break spatial symmetry and favor one or two 
dimensions. As illustrated in the right panel of Fig. [3l 
for example, multiple 'pockets' can coalesce into locally 
or even globally (e.g., stripe-like) connected structures, in 
contrast with the Wigner crystal state where each pocket, 
corresponding to one electron, is fully localized. 

To conclude, we have determined the true HF ground 
state for finite electron gas. Combining numerical and 
analytic results, we have described the origin and char- 
acteristics of the broken symmetry state at high density, 
and the novel pairing mechanism that drives it. 
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